#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>

#include <algorithm>
#include <vector>

using std::vector;
using std::sort;

typedef long int N;
typedef long double F;
typedef F TBarCubaturPoint[4]; // a1 a2 a3 w
typedef F TCubaturPoint[3]; // x y w

// global variables defined from console parameters
const TBarCubaturPoint* barformula;
TCubaturPoint formula[128];
int count;
F EPS;
N MAX_L;
F B[6];
F &X1 = B[0], &Y1 = B[1], &X2 = B[2], &Y2 = B[3], &X3 = B[4], &Y3 = B[5];
F Det;
F Square;

void Transform(const TBarCubaturPoint& a, TCubaturPoint& p, F B[] = B) {
  p[0] = a[0]*B[0] + a[1]*B[2] + a[2]*B[4];
  p[1] = a[0]*B[1] + a[1]*B[3] + a[2]*B[5];
  p[2] = a[3];
}

N HalfPow2(N n) {
  if (n & (n - 1)) {
    while (int k = n & (n - 1))
      n = k;
  } else
    n = n / 2;
  return n;
}

struct TSummator {
  vector<F> table;
  TSummator (F f = 0.0) : table(1, f) {}
  void operator += (F f) {
    if (table.size() >= 1024 * 1024)
      calc();
    table.push_back(f);
  }
  F calc() {
    sort(table.begin(), table.end());
    while (table.size() > 1) {
      N p = HalfPow2(table.size());
      for (int i = 0; i < p; ++i)
        table[i] = (2*i   < table.size()? table[2*i]  : 0.0) + 
                   (2*i+1 < table.size()? table[2*i+1]: 0.0);
      table.resize(p);
    }
    table.assign(1, table[0]);
    return table[0];
  }
  operator F () { return calc(); }
} sum;

F Power(F x, N n) {
  F res = 1.0;
  while (n) {
    if (n & 1)
      res *= x;
    x *= x;
    n = n >> 1;
  }
  return res;
}

F CubaturCalc(int r, int s) {
  TSummator sum = 0;
  for (size_t i = 0; i < count; ++i)
    sum += formula[i][2] * Power(formula[i][0], r) * Power(formula[i][1], s);
  return sum;
}

F CubaturCalc(int r, int s, F triangle[]) {
  for (size_t i = 0; i < count; ++i)
    Transform(barformula[i], formula[i], triangle);
  F sum = 0;
  for (size_t i = 0; i < count; ++i)
    sum += formula[i][2] * Power(formula[i][0], r) * Power(formula[i][1], s);
  return sum;
}

N GCD(N a, N b) {
  while (a && b) {
    if (a < b)
      b %= a;
    else 
      a %= b;
  }
  return a? a: b;
}

F IntegrateBar(int i, int j, int k) {
  F res = 1.0;
  int c = 1, z = 2, i0 = 1, j0 = 1, k0 = 1, z0 = 2;
  while (i0 <= i || j0 <= j || k0 <= k) {
    int zdelta = (i0 <= i) + (j0 <= j) + (k0 <= k);
    int m = (i0 <= i? i0++: 1) * (j0 <= j? j0++: 1) * (k0 <= k? k0++: 1);
    if (c > (1 << 20) / m || z > (1 << 20) / z0) {
      res *= (F)c / (F)z;
      c = z = 1;
    }
    while (zdelta--)
      z *= ++z0;
    c *= m;
    N gcd = GCD(c, z);
    c /= gcd;
    z /= gcd;
  }
  return res * (F)c / (F)z;
}

F IntegrateTriangle(int r, int s, F B[]) {
  F &X1 = B[0], &Y1 = B[1], &X2 = B[2], &Y2 = B[3], &X3 = B[4], &Y3 = B[5];
  TSummator sum = 0;
  for (int i = 0; i <= r; ++i)
    for (int j = 0; j <= r - i; ++j)
      for (int k = 0; k <= s; ++k)
        for (int p = 0; p <= s - k; ++p) {
          F pwr = Power(X1, i) * Power(X2, j) * Power(X3, r - i - j) *
                  Power(Y1, k) * Power(Y2, p) * Power(Y3, s - k - p);
          sum += 2 * pwr * IntegrateBar(i + k, j + p, r - i - j + s - k - p);
        }
  return sum;
}

F IntegrateTriangle(int r, int s, int step, F B[]) {
  F &X1 = B[0], &Y1 = B[1], &X2 = B[2], &Y2 = B[3], &X3 = B[4], &Y3 = B[5];
  sum = 0;
  if (step) {
    F triangles[4][6];
    triangles[0] = {X1, Y1, (X1+X2), (Y1+Y2)/2.0, (X1+X3)/2.0, (Y1+Y3)/2.0};
    triangles[1] = {X2, Y2, (X2+X3), (Y2+Y3)/2.0, (X2+X1)/2.0, (Y2+Y1)/2.0};
    triangles[2] = {X3, Y3, (X3+X1), (Y3+Y1)/2.0, (X3+X2)/2.0, (Y3+Y2)/2.0};
    triangles[3] = {(X1+X2)/2.0, (Y1+Y2)/2.0, (X2+X3)/2.0, (Y2+Y3)/2.0, (X1+X3)/2.0, (Y1+Y3)/2.0};
    return (IntegrateTriangle(r, s, step - 1, triangles[0]) +
            IntegrateTriangle(r, s, step - 1, triangles[1]) +
            IntegrateTriangle(r, s, step - 1, triangles[2]) +
            IntegrateTriangle(r, s, step - 1, triangles[3]))/4.0;
  } else
    //return IntegrateTriangle(r, s, B);
    return CubaturCalc(r, s, B);
}

#define ARRAY_SIZE(a) sizeof(a)/sizeof(*(a))

#define SYMM6(w, a, b, c) {a, b, c, w}, {a, c, b, w}, {c, a, b, w}, \
                          {c, b, a, w}, {b, c, a, w}, {b, a, c, w}
#define SYMM3(w, a, b)    {a, a, b, w}, {a, b, a, w}, {b, a, a, w}
#define SYMM1(w, a)       {a, a, a, w}
#define FORMULA(index) TBarCubaturPoint formula##index[]
#define _FORMULA(index) formulas[index] = formula##index, counts[index] = ARRAY_SIZE(formula##index);
const TBarCubaturPoint* formulas[10];
int counts[10];

  FORMULA( 0) = {SYMM3( 0.333333333333333, 0.166666666666667, 0.666666666666666)};
  FORMULA( 1) = {SYMM3( 0.333333333333333, 0.500000000000000, 0.000000000000000)};
  FORMULA( 2) = {SYMM1(-0.562500000000000, 0.333333333333333),
                 SYMM3( 0.520833333333333, 0.200000000000000, 0.600000000000000)};
  FORMULA( 3) = {SYMM6( 0.166666666666667, 0.659027622374092, 0.231933368553031, 0.109039009072877)};
  FORMULA( 4) = {SYMM3( 0.109951743655322, 0.091576213509771, 0.816847572980458),
                 SYMM3( 0.223381589678011, 0.445948490915965, 0.108103018168070)};
  FORMULA( 5) = {SYMM1( 0.375000000000000, 0.333333333333333),
                 SYMM6( 0.104166666666667, 0.736712498969535, 0.237932366475534, 0.025355134554931)};
  FORMULA( 6) = {SYMM3( 0.205950504760887, 0.437525248383384, 0.124949503233232),
                 SYMM6( 0.063691414286223, 0.797112651860071, 0.165409927389841, 0.037477420750088)};
  FORMULA( 7) = {SYMM1( 0.225000000000000, 0.333333333333333),
                 SYMM3( 0.125939180544827, 0.101286507323456, 0.797426985353088),
                 SYMM3( 0.132394152788506, 0.470142064105115, 0.059715871789770)};
  FORMULA( 8) = {SYMM3( 0.050844906370207, 0.063089014491502, 0.873821971016996),
                 SYMM3( 0.116786275726379, 0.249286745170910, 0.501426509658180),
                 SYMM6( 0.082851075618374, 0.636502499121399, 0.310352451033785, 0.053145049844816)};
  FORMULA( 9) = {SYMM1(-0.149570044467670, 0.333333333333333),
                 SYMM3( 0.175615257433204, 0.260345966079038, 0.479308067841924),
                 SYMM3( 0.053347235608839, 0.065130102902216, 0.869739794195568),
                 SYMM6( 0.077113760890257, 0.638444188569809, 0.312865496004875, 0.048690315425316)};

int rs_rows;
int r_max;
F eps_max = 0;
vector<vector<F>> rs_table;
vector<F> eps_table[2];

void build_rs_table() {
  bool zero = true;
  rs_table.resize(0);
  for (int n = 0; zero; ++n) {
    rs_table.resize(rs_table.size() + 1);
    for (int r = 0; r <= n; ++r) {
      F I = IntegrateTriangle(r, n - r, B);
      F S = CubaturCalc(r, n - r);
      rs_table.back().push_back(fabs((I-S)/I));
      zero = zero && fabs(rs_table.back().back()) < EPS;
      //printf("%lg(I=%lg,S=%lg),", (double)rs_table.back(),I,S);
      if (eps_max < fabs(rs_table.back().back())) {
        eps_max = fabs(rs_table.back().back());
        r_max = r;
      }
    }
    rs_rows = n + 1;
  }
      //printf("\n");
}

void build_eps_table() {
  eps_table[0].resize(MAX_L);
  eps_table[1].resize(MAX_L);
  F prev = rs_table[rs_rows - 1][r_max];
  //F S = CubaturCalc(r_max, rs_rows - r_max);
  F I = IntegrateTriangle(r_max, rs_rows - r_max, B);
  for (int l = 0; l < MAX_L; ++l) {
    F S = IntegrateTriangle(r_max, rs_rows - r_max, l + 1, B);
    eps_table[0][l] = fabs((I-S)/I);
    eps_table[1][l] = 1 + 1/((F)l*l*l*l);
    prev = eps_table[0][l];
  }
}

void print_rs_table() {
  printf("\\begin{tabular}{|");
  for (int i = 0; i <= rs_rows; ++i)
    printf("c|");
  printf("}\n\\hline\n$r\\backslash s$");
  for (int i = 0; i < rs_rows; ++i)
    printf(" & $%d$", i);
  printf(" \\\\ \\hline\n");
  int n = 0;
  for (int r = 0; r < rs_rows; ++r) {
    printf("$%d$", r);
    for (int s = 0; s < rs_rows; ++s)
      if (s < rs_rows - r)
        printf(" & $%lg$", rs_table[r+s][r]);
      else
        printf(" & ... ");
    printf(" \\\\ \\hline\n");
  }
  printf("\\end{tabular}\n");
}

void print_eps_table() {
  printf("\\begin{tabular}{|");
  for (int i = 0; i <= MAX_L; ++i)
    printf("c|");
  printf("}\n\\hline\n$l$");
  for (int i = 0; i < MAX_L; ++i)
    printf(" & $%d$", i + 1);
  printf(" \\\\ \\hline\n");
  int n = 0;
  printf("$\\epsilon_{%d, %d}^l$", r_max, rs_rows - r_max);
  for (int i = 0; i < MAX_L; ++i)
    printf(" & %lg", eps_table[0][i]);
  printf(" \\\\ \\hline\n");
  printf("$\\epsilon_{%d, %d}^l \\slash \\epsilon_{%d, %d}^{l-1}$", 
         r_max, rs_rows - r_max, r_max, rs_rows - r_max);
  for (int i = 0; i < MAX_L; ++i)
    printf(" & %lg", eps_table[1][i]);
  printf(" \\\\ \\hline\n");
  printf("\n\\end{tabular}\n");
}

F TaskCalc(F triangle[]) {
  //printf("TRIANGE: %lg %lg %lg %lg %lg %lg\n", triangle[0], triangle[1], triangle[2], triangle[3], triangle[4], triangle[5]);
  for (size_t i = 0; i < count; ++i) {
    //printf("Bar: %lg %lg %lg %lg\n", barformula[i][0], barformula[i][1], barformula[i][2], barformula[i][3]);
    Transform(barformula[i], formula[i], triangle);
  }
  F sum = 0;
  for (size_t i = 0; i < count; ++i) {
    //printf("Calc: %lg %lg %lg = %lg *%lg\n", formula[i][0], formula[i][1], formula[i][2], 
    //               formula[i][2] * sqrt(formula[i][0]) * exp(-formula[i][0]*formula[i][1]),
    //               -formula[i][0]*formula[i][1]);
    sum += formula[i][2] * sqrt(formula[i][0]) * exp(-formula[i][0]*formula[i][1]);
  }
  return sum;
}

F Task() {
  F triangle[6];
  F S = 0;
  const F dist = 8.8;
  const int x_steps = 1000;
  const int y_steps = 32;
  for (int i = 0; i < x_steps; ++i) {
      for (int j = 0; j <= y_steps; ++j) {
          F x[4] = {(F)j/(F)y_steps, (F)j/(F)y_steps, (F)(j+1)/(F)y_steps, (F)(j+1)/(F)y_steps};
          F y[4] = {i*dist, (i+1)*dist, (i+1)*dist, i*dist};
          triangle = {x[0], y[0], x[1], y[1], x[2], y[2]};
          S += (dist*TaskCalc(triangle))/((F)2*y_steps);
          triangle = {x[0], y[0], x[2], y[2], x[3], y[3]};
          S += (dist*TaskCalc(triangle))/((F)2*y_steps);
      }
  }
  return S;
}

void print_tables() {
  printf("\\documentclass[draft, a4paper, 12pt]{article}\n");
  printf("\\begin{document}\n");
  print_rs_table();
  printf(" \\\\\n\n\n");
  print_eps_table();
  printf(" \\\\\n\n\n");
  printf("Task: %lg\n", (double)Task());
  printf("\\end{document}\n");
}

int main(int argc, char* argv[]) {
  if (argc != 9 && argc != 10) {
    if (argc != 1)
      printf("Error: incorrect count of arguments\n");
    printf("Usage: prog <formula> <eps> <x1> <y1> <x2> <y2> <x3> <y3> [MAX_L]\n"
           "  formula -- integer number in table\n"
           "  eps     -- float point precision for comparation\n"
           "  (x1 y1) (x2 y2) (x3 y3) -- float point coordinates of triangle\n"
           "  MAX_L   -- [default: 8] rows in table 2\n");
    return 0;
  }

  _FORMULA( 0); _FORMULA( 1); _FORMULA( 2); _FORMULA( 3); _FORMULA( 4); _FORMULA( 5);
  _FORMULA( 6); _FORMULA( 7); _FORMULA( 8); _FORMULA( 9);

  int index = atoi(argv[1]) - 1;
  EPS = atof(argv[2]);
  for (int i = 0; i < 6; ++i)
    B[i] = atof(argv[3 + i]);
  MAX_L = argv[9]? atoi(argv[9]): 8;

  count = counts[index];
  barformula = formulas[index];
  for (int i = 0; i < count; ++i)
    Transform(barformula[i], formula[i]);

  F Det = B[0]*B[3] + B[2]*B[5] + B[1]*B[4] 
        - B[3]*B[4] - B[1]*B[2] - B[0]*B[5];
  F Square = Det / 2.0;

  build_rs_table();
  build_eps_table();
  print_tables();

  /*printf("count = %d\n", count);
  sum = 0;
  sum += 1;
  sum += 2;
  sum += 3;
  sum += 4;
  sum += 5;
  sum += 6;
  sum += 7;
  sum += 8;
  sum += -6;
  printf("hp2 %d %d\n", 7, HalfPow2(7));
  printf("hp2 %d %d\n", 6, HalfPow2(6));
  printf("hp2 %d %d\n", 5, HalfPow2(5));
  printf("hp2 %d %d\n", 4, HalfPow2(4));
  printf("sum = %lf\n", (double)sum);
  printf("2.0^3 = %lf\n", (double)Power(2.0, 3));
  printf("2.01^5 = %lf\n", (double)Power(2.01, 5));
  printf("3^4 = %lf\n", (double)Power(3.0001, 4));
  printf("CubaturCalc: %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg\n", 
                     (double)Square*CubaturCalc(0, 0),
                     (double)Square*CubaturCalc(1, 0),
                     (double)Square*CubaturCalc(0, 1),
                     (double)Square*CubaturCalc(2, 0),
                     (double)Square*CubaturCalc(1, 1),
                     (double)Square*CubaturCalc(0, 2),
                     (double)Square*CubaturCalc(3, 0),
                     (double)Square*CubaturCalc(2, 1),
                     (double)Square*CubaturCalc(1, 2),
                     (double)Square*CubaturCalc(0, 3)
                );
  printf("IntegrateBar: %lg\n", 
           //(double)IntegrateBar(0, 0, 0)//,
           //(double)IntegrateBar(1, 0, 0)//,
           (double)IntegrateBar(2, 3, 3)//,
           //(double)IntegrateBar(1, 1, 0)//,
           //(double)IntegrateBar(1, 1, 1)//,
           //(double)IntegrateBar(2, 1, 1)
        );
  printf("Integrate: %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg\n", 
           (double)Square*IntegrateTriangle(0, 0, B),
           (double)Square*IntegrateTriangle(1, 0, B),
           (double)Square*IntegrateTriangle(0, 1, B),
           (double)Square*IntegrateTriangle(2, 0, B),
           (double)Square*IntegrateTriangle(1, 1, B),
           (double)Square*IntegrateTriangle(0, 2, B),
           (double)Square*IntegrateTriangle(3, 0, B),
           (double)Square*IntegrateTriangle(2, 1, B),
           (double)Square*IntegrateTriangle(1, 2, B),
           (double)Square*IntegrateTriangle(0, 3, B)
        );*/
}
